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Abstract. A mean-field sandpile model that exhibits self-organized criticality (SOC) despite violation of 
the grain-transfer conservation law during avalanches is proposed. The sandpile consists of A'^ agents and 
possesses background activity with intensity rj G [0, 1]. Background activity is characterized by transitions 
between two stable agent states. Analysis employing theories of branching processes and fixed points reveals 
a transition from sub-critical to SOC phase that is determined by rjN. The model is used to explain the 
school size distribution of free-swimming tuna as a result of population depletion. 
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1 Introduction 

Prototypical SOC sandpile models have attracted several 
researchers due to their inherent simplicity yet non-trivial 
behavior. Dynamical mean-field theory has been proven to 
offer the most effective elucidation of how SOC works 1 . 
One of the key ingredients found to bring about SOC 
in sandpiles is grain-transfer conservation. Tsuchiya and 
Katori offered rigorous proof that violation of the con- 
servation law frustrates the criticality of Abelian sand- 
piles A mean-field treatment of the Manna sandpile 
known as the self-organized branching process (SOBP) 
has been proposed, further demonstrating that noncon- 
servation disrupts criticality [3,. Breaking of SOC due 
to nonconservation has also been later demonstrated nu- 
merically and analytically in the Olami-Feder-Christensen 
earthquake model [4j. The forest-fire model (FFM), which 
does not have conservation laws, was also proposed as a 
model for SOC [5]. However, it was ultimately proven, via 
analysis of Lyapunov exponents [6j and later by means of 
renormalization group approach [7], that FFM does not 
exhibit SOC. 

Signatures of scale invariance and criticality in biology, 
ecology, and a wide range of other animate complex sys- 
tems have been linked to the principles behind power-law 
generating SOC models such as the sandpile model [8^. 
Practically, biological systems are not expected to be con- 
servative either because they are constantly interacting 
with the environment and are far from equilibrium. These 
systems are therefore not supposed to display SOC. 

Here, I propose a mean-field sandpile model that dis- 
plays criticality despite violation of grain-transfer conser- 
vation. This model takes off from the conceptual difficulty 
with typical sandpile models in defining fluctuations be- 



tween stable states because passive sites only switch states 
when an avalanche crosses them or when a grain is sprin- 
kled on them. Thus, sites in this model are turned into 
active agents. Demonstrating criticality in self-organizing 
systems that are nonconservative enhances the synthetical 
capacity of SOC theory in the arena of biocomplexity. 

Section [2] describes the nonconservative model. Sec- 
tion [3] presents and discusses the results of analysis and 
simulations of the model. Section |4] elaborates salient as- 
sumptions and possible areas for further analysis and ex- 
tension. Lastly, Section [5] summarizes and concludes the 
paper. 



2 Model 

2.1 Framework 

Inspired by biocomplexity, a sand grain assumes the role of 
a stimulus — for instance an environmental trigger — that 
is stored, integrated, and transferred from one agent to 
another. The build-up of stored stimuli induces agent re- 
sponse, which activates the transfer of stimulus through 
neighboring agents. Stimulus propagation corresponds to 
avalanches. 

In accordance with Manna sandpile rules and studies 
of excitable media [3], each agent in the system may be 
in any three states {refractory, quiescent, excited} mapped 
to {z} = {0,1,2}, as follows: 

• quiescent z — 1 

• excited i-^ z = 2 

• refractory i-^ z — 
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By storage and integration of stimulus, refractory agents 
may turn quiescent and quiescent agents may become ex- 
cited. If threshold zth — 1, then an excited agent's stim- 
ulus level z > zth, so that it fires a response in order to 
relax back to the stable states z = 0, 1. When this hap- 
pens, an avalanche ensues. "Avalanche" is here interpreted 
as clustering in the sense that the agents excited during 
an avalanche are all "behaviorally matched" to the initial 
agent that triggers the avalanche. 

This model differs from typical sandpile models be- 
cause it allows for the possibility that, apart from excita- 
tions triggered along the path of the avalanche, a refrac- 
tory agent can spontaneously become quiescent and vice 
versa, albeit at a much slower pace than the avalanche it- 
self. This independent but slow process is what is referred 
hereafter as background activity. 

The system is a population of iV = 2"+^ — 1 excitable 
agents, where n corresponds to the maximum number of 
generations of excitation allowed during a single avalanche, 
so that n is a boundary condition. The model has no 
explicit dependence on any spatial lattice configurations, 
which makes analysis straightforward due to the absence 
of nearest-neighbor correlations and is a more accurate 
representation of socially interacting systems made up of 
entities that are constantly in motion such as animal groups. 
The mean-field nature of the model also does not require 
a physical boundary. Fish schools, for instance, traverse a 
seemingly limitless oceanic space. 



2.2 Mean-Field Transition Rules 

The population dynamically reconfigures the landscape of 
agent states when there is no avalanche. A quiescent agent 
turns refractory with probability A, whereas a refractory 
agent turns quiescent with probability 77. These transi- 
tions between refractory and quiescent states (stable agent 
states) is the background activity. 

The time-dependent density q{t) of quiescent agents 
also corresponds to the probability that an external stim- 
ulus injected at time t into a randomly selected agent initi- 
ates excitation, because only quiescent agents can become 
excited at any given time. If at time t an excited agent in- 
deed emerges, time is frozen and avalanche ensues. Freez- 
ing time follows from the assumption that background ac- 
tivity takes place much slower than an avalanche — known 
as "timescale separation," which is a typical assumption 
in sandpile models jl]. 

During an avalanche, an excited agent stimulates at 
most two other agents. With probability a the excited 
agent turns refractory by transferring 2 units of stimuli 
to two randomly chosen neighbors. With probability (3 
the excited agent turns quiescent by dissipating only 1 
unit of stimulus to one random neighbor. Incorporating 
a probability e = I — a — (3 that an excited agent turns 
refractory by inwardly absorbing 2 units of stimuli without 
subsequently stimulating other agents, the transfer rule 
becomes nonconservative. Conservation law is therefore 
violated for e > 0. 



Table 1. Mean- field transition rules of the nonconservative 
branching model with corresponding transition probabilities. 
The central digits for the avalanche rules correspond to the ex- 
cited agent and the flanking digits are the two random neigh- 
bors, in no particular order. The sum of transition probabilities 
for the avalanche and nonconservative rules is q, which is, self- 
consistently, the probability that the avalanche had initiated. 
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2^0 




eg = (1 - Q - I3)q 



The nonconservative transfer rule repeats until all ex- 
cited agents are depleted, after which the frozen time takes 
off again and the population builds up for the next avalanche. 
There is a limit n to the number of generations of transfers 
during an avalanche. At the n-th generation, any remain- 
ing excited agents mandatorily absorbs stimuli without 
further stimulation of other agents. All the aforementioned 
processes are summarized in table [1] 



2.3 Dynamics and Branching Process 

Given the transition rules, the density g, assumed to be 
a continuous dynamical variable, satisfies the stochastic 
dynamical equation 

'^^={l-q)7^-q\ + A{q;a,P)+miN. (1) 

The noise term S^/N arises from the stochastic transition 
rules, accounting for fluctuations around mean values as- 
sumed to hold in the mean-field calculations. It appropri- 
ately vanishes in the large-iV limit. The term A{q\ a, (3) 
represents the change in q{t) due to redistribution of stim- 
ulus by an avalanche at time t. Treating an avalanche as a 
branching process, and following closely the analysis in [3], 
A{q] a, P) satisfies 



l-il-e)q 



1 



1 



1 



-2ct' 



, (2) 



where a is the branching parameter derived from the fol- 
lowing definition 

(T = ^fc7rfc, (3) 

k 

wherein from table [U the branching probability tt^ that 
an excited agent subsequently stimulates k other agents is 



TTfc = aqSk,2 + PqSk.i + [!-(!- e)q] Sk^o, 



(4) 



with dij being the Kronecker delta. The first and sec- 
ond terms in eq. ^ are the sum of transition probabili- 
ties of the first three avalanche rules and of the last three 
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avalanche rules listed in table[l] respectively. The last term 
is the total probability coming from the nonconservation 
rule and the case wherein a subsequently stimulated agent 
is refractory, which both result to fc = 0. Substitution of 
eq. dH) into eq. ^ gives 



cr = (2q + I3)q, 



(5) 



which depends on a and /3, and is proportional to den- 
sity q. In this model it is assumed that a and /3 remain 
fixed throughout so that the only dynamic variable is q. 
Consequently, 

P)t (6) 



dt 



(2a 



dt 



A branching process is sub-critical when cr < 1, and 
consequently avalanches have sizes always smaller than a 
finite cutoff size. On the other hand, a branching process 
is supra-critical when cr > 1, and consequently avalanches 
with sizes as large as the system itself are formed almost 
with certainty. Hence, ct = 1 is a critical value that results 
to a critical branching process [lOj . The avalanche size 
distribution emerging from a critical branching process 
is expected to be power law. Since a only varies when q 
changes, the q determines whether the branching process 
is sub-critical, critical, or supra-critical. 

To keep the model minimal, the number of parameters 
in eq. ([1]) is reduced by introducing the relation 



A 



a/q-1 
2a + 13-1, 



(7) 



which effectively couples the background activity to the 
avalanche (branching process). Eq. ([7]) is a phenomeno- 
logical assumption which essentially guarantees that q(t) 
in the steady state approaches a value that makes a — 1. 
This is an indication that the model self-organizes to its 
critical state. In order to analytically prove that the model 
indeed exhibits SOC, the critical state must be an attrac- 
tive fixed point in phase space [71. The phase portrait is 
simply a plot of q := dq/dt versus q. Since a (x q, then the 
phase portrait can equivalently be visualized by plotting 
a := dcr/dt, defined in eq. (O, versus a. 



2.4 Avalanche Size Distribution 

An avalanche resulting from the rules of the model is illus- 
trated as a branching tree in figure [TJ The initial excited 
agent is the topmost node, giving rise to two subsequently 
excited nodes with probability aq. One of these excited 
nodes in turn generates one excited node with probabil- 
ity Pq, while the other generates two excited nodes with 
probability aq. One of such excited nodes do not gener- 
ate a subsequently excited node with probability eq due 
to nonconservation. Setting n = 3, all excited nodes at 
three levels below the topmost node do not further gen- 
erate excited nodes. In the example shown, the avalanche 
size s = 9. 




Fig. 1. Avalanche for n — 3 shown as a branching tree. Shaded 
circles correspond to excited agents. The avalanche size s — 9 
corresponds to the total number of shaded circles. 



The branching probability iTk defined in eq. ([4]) essen- 
tially describes the likelihood that an excited agent sub- 
sequently generates k g {0,1,2} excited agents in the 
succeeding generation. Using tt^ the avalanche size dis- 
tribution P{s) is calculated via a generating functional 
formalism. A generating function 5'm(w) for P(s) after m 
generations is defined as 



5m (w) 



(8) 



Incidentally, 'Sm is also the m-th iterate of the generating 
function := 5, defined as 



E 

s=l 



(9) 



Also by definition, ^m+i — d{dm), such that from eq. ^ 
one obtains 



l?m+l(w) 

which simplifies to 



: UJ 



{aqSl+(3qS,n + [l-il~e)q]} , (10) 



following from the definition of tt^ in eq. (|4|). For large 
enough m, the theory of branching processes asserts a self- 
consistency relation wherein ^m+i — i?m, so that solving 
for in eq. pO|) yields 



i?m(w) 



\ — buj — Vl — 26w -|- aw^ 
2aquj ' 



(11) 



where a — 0^q^ — 4ag[l — (1 — e)?], and h — f3q. Binomial 
expansion of eq. (jlip around its singularity uj = 0, one 
obtains a power series similar to eq. The coefficients 
of this expansion correspond to P{s) from eq. ([8|). In a 
more compact form the solution can be expressed as a 
recurrence relation 



P{.s) = 



1 



1 



[(2s - l)bP{s - 1) - (s - 2)aP{s - 2)] , 



(12) 
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Fig. 2. Avalanche size distribution of a nonconservative system 
with e = 0.2 and rjN « 8192. A power-law with exponent 
Ts =3/2 is drawn as guide to the eye. The exponential fat tail 
is an artifact of the finiteness of the system. The inset graph is 
the distribution of avalanche lifetimes, which is the number of 
excitation generations before the avalanche dies out. A power 
law with exponent rx = 2 is drawn as guide. 



Eq. (jl2[) may be shown, by graphical inspection, to have 
the asymptotic behavior P{s) ~ s~'^/^ for s '3> 1 if a = 
26 — 1 ; a condition that can equivalently be expressed as 



Q{q) = a + 25- 1 

= /32g2 _ 4aq[l - (1 - e)q] -2f3q+l 
= 0. 



(13) 



Numerical simulation of the model also confirms the asymp- 
totic behavior of eq. p2| , as shown in figure [2l for degree 
of nonconservation e — 0.2 and rjN « 8192. 

The function Q{q) is parabolic in terms of q and has a 
unique root at qc = (2a+ , which makes a = 1. Hence, 
eq. (|13p is a criticality condition. Furthermore, since a cx 
q, then Q{q) can alternatively be expressed as Q(cr) which 
has a root at a = 1. 



3 Results and Discussion 

The stationary behavior of the model is examined using 
a well-established geometrical theory of fixed points [llj . 
Nonlinear differential equations such as eq. ([6]) may be an- 
alyzed graphically in terms of vector fields. In this frame- 
work & is interpreted as a "velocity vector" at each possi- 
ble a value. Plotting a versus a is the phase portrait of the 
model. The fixed point a* is the value of ct at which a — 0. 
The trajectory of the vector around the neighborhood of 
this fixed point is directed to the right where tr > 0, and 
to the left where a < 0. This means that if & is increasing 
around the neighborhood of a*, then the fixed point is 
repulsive. On the other hand, if a is decreasing, then the 
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Fig. 3. Phase portrait showing cr versus cr, of a nonconservative 
system with e = 0.25 and r\ = 0.03125 for different n: 4 (chain 
curve), 8 (dotted curve), 16 (dashed curve), and oo (broken 
line). The root of a corresponds to the fixed point value a* 
of the branching parameter. The solid parabolic curve is the 
function Qicr) with root at cr = 1. For n > 16, a* — 1, which 
indicates that the nonconservative system evolves towards its 
critical state via self-organization. 



fixed point is attractive. Linearization of eq. ([6]) actually 
yields 

,• 9 da ^ 
hm —— < 0. 

o-^CT* oa at 

Shown in figure [3] are the phase portraits of the model 
for different n, which are monotonically decreasing. Thus, 
the fixed points of eq. ([6|) are indeed dynamically attrac- 
tive. Therefore, apart from fluctuations, the model spon- 
taneously approaches the state defined by its fixed point — 
the very idea of self-organization. 

Also illustrated in figure [3] is the function Q{(j) with a 
root at cr = 1. For n = A and n = 8, obviously cr* < 1, im- 
plying that the system self-organizes towards a sub-critical 
state. However, from n = 16 up to n —> cx), | cr* — 1 | ^ 1. 
A nonconservative system achieves criticality even if its 
size is finite. The rapid approach of cr* towards 1 from 
n = 8 to n = 16 and the decelerating change in cr* from 
n = 16 to n — > oo indicates a phase transition. 

Indeed, as shown in figure 31 a phase transition of 
the model exists. The data points are derived from nu- 
merical fixed-point analysis of eq. ([1]), assuming that the 
noise term (^/N is very small to be significant. For a wide 
range of parameter values: e g (0, 1/2) and 77 g (0, 1], a 
transition surfaces out from the relation between rjN and 
(g*— 1/2) / ((7c— 1 /2). The resultant plot is fitted by a logistic 
curve 

q* - 1/2 



qc - 1/2 



1 



1 



TjN 



(14) 



where 9 

,2 



9.31 ± 0.27 and 7 

-5 Tl2 



0.68 ± 0.01 (goodness of 
fit: xVDoF = 4 X 10-^ B? = 0.99986, no weighting). 
Eq. (fH)) suitably describes two limiting cases. The first 
one is ?7 ^ 0, which results to q* 1/2. This limiting 
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Fig. 4. Phase transition in the nonconservative branching 
modeL Data points result from numerical analysis of the fixed 
points q* of Equation [T] neglecting fluctuations for var- 

ious e € (0,0.5), Tj G (0,1], and system size TV. The curve 
is a logistic fit of the form y{x) = 1 — [1 + {x/6)'^\~^ where 
y = {q* - l/2)/(gc - 1/2) and x = t)N . 



case is equivalent to having no background activity, as 
in the SOBP model by Lauritsen, Zapperi & Stanley [3 . 
They have similarly found that the steady-state (fixed- 
point) value is q* = 1/2, irrespective of any degree of 
nonconservation e. Consequently, their sandpile model is 
al"ways sub-critical whenever e > 0. 

The second limiting case is nN ^ 9, which leads to 
q* — » gc- At this limit, the nonconservative system is al- 
ways critical. A closer look at figured] reveals that the data 
points lie considerably close to 1 starting at rjN ~ lO**. 
This means that if the population is very large, say N = 
10^, then for a wide range of r/ £ (10~^, 1] the nonconser- 
vative system achieves SOC. Simulation of the model for 
various degrees of nonconservation e yields an avalanche 
size distribution of the form plotted in figure [21 

In the context of ecological systems, large N may easily 
be satisfied by clumps of microorganisms such as cellular 
aggregates in culture [12j . However, in macroscopic sys- 
tems such as animal groups, this may not be the case be- 
cause of several ecological factors, such as the presence of 
predators and destruction of habitats, which curtail the 
proliferation of a certain species. A good example is a 
population of tuna fishes studied by Bonabeau et al. p[3|. 
The school size distribution P{s) found from this study 
is plotted in fig. [5l fitted by the nonconservative model 
with e = 0.15 and r]N = 409.5. This rjN value maps to 
a sub-critical phase on the phase diagram shown in fig- 
ure [H explaining the prominent exponential tail in P{s). 
Also shown is a curve from an equivalent SOBP model 3 , 
exhibiting less agreement with the data. 

A school size distribution having an exponential tail 
implies that schools cannot be expected to be larger than 
a finite cutoff size. However, large schools or groups of 
animals are required for sustainability and wildlife preser- 
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Fig. 5. School size distribution of free-swimming tuna (shaded 
circles), and of a nonconservative system with e — 0.15 and 
rjN = 409.5 taken after 2^^ time steps (solid curve). Also shown 
is a size distribution resulting from an equivalent SOBP model 
(dashed curve) having the same value of Qc = 0.952 as the 
nonconservative model. Inset graph reveals the evolution of 
q{t). Its steady-state value q* is less than qc (dashed line), 
implying that the nonconservative system is sub-critical. 

vation as it is believed that group size affects reproductive 
success of animals, most especially endangered species. 

Lastly, it is also interesting to note that for rjN ~ 9, the 
q* value is highly uncertain to small variations in either rj 
or N. This is the regime describing the abrupt rise from 
to 1 in figure m 



4 Recommendations 

A key assumption in the model rests upon having the 
branching parameter defined in eq. ^ sufficiently close 
to its critical value of 1. With a and /3 fixed, cr = 1 is ob- 
tained solely by the incorporation of a background activ- 
ity. This background activity is coupled to the avalanche 
via a matching condition [eq. ([7])] to guarantee that q = 
{2a + f3)~^ at the steady state. As a modification, one 
could assume instead that both a and (3 also change with 
t, perhaps due to some local feedback mechanism that ad- 
justs the branching probability depending on the history 
of an agent's activity. Hence in general eq. ([6]) may be 
written as 

do- /„5q; 9/3\ ,^ ^Bq 

d^n'^^^)'^^^'"^^^^- 

The implementation of the above modification lies in defin- 
ing the feedback mechanism that gives da/dt and df3/dt. 
Lastly, the model can be extended to the case wherein 
{z} = {0, 1, 2, ... , zth, 2th + 1} and zth > 1- Preliminary 
research is underway to evaluate the performance of the 
model with these modifications. 
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5 Summary and Conclusion 

In summary, I have proposed a mean-field sandpile model 
coupled with a baekground activity that is characterized 
by local fluctuations between stable agent states. Even 
in the presence of violation of the grain-transfer rule, the 
sandpile can self-organize to a critical state. This result 
addresses the long-standing issue of whether SOC is pos- 
sible in nonconservative sandpiles. The model is applied to 
explaining the truncated school size distribution observed 
for free-swimming tuna and provides insight into the ef- 
fects of population depletion on the aggregation capacity 
of clustering animals. Lastly, some recommendations for 
extending the model have been rendered. 

The author wishes to acknowledge the OVCRD in UP Diliman 
for Grant No. 050501 DNSE, and the PCASTRD-DOST for 
funding this research. 
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